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INTEGRATION  OF  STEREO  VISION  AND  OPTICAL  FLOW 
USING  ENERGY  MINIMIZATION  APPROACH 

Nasser  M.  Nasrabadi 

Computer  Vision  Research  Group 
Electrical  Engineering  Department 
Worcester  Polytechnic  Institute 
Worcester,  Massachusetts  01609 
Tel:  (716)  636-2427 

ABSTRACT 

A  cooperative  motion-stereo  method  is  proposed  where  image  intensity 
(brighmess)  and  optical  flow  information  are  integrated  into  a  single  stereo 
technique  by  modeling  the  input  data  as  coupled  Markov  Random  Fields 
(MRF).  The  Bayesian  probabilistic  estimation  method  and  the  MRF-Gibbs 
equivalence  theory  are  used  to  integrate  the  optical  flow  and  the  gray  level 
intensity  information  to  obtain  an  energy  function  which  will  explicitly 
represent  the  depth  discontinuity  and  occlusion  constraints  on  the  solution. 
This  energy  function  involves  the  similanty  in  intensity  (or  edge  orientation) 
and  the  optical  flow  between  corresponding  sites  of  the  left  and  right  images 
as  well  as  the  smoothness  constraint  on  the  disparity  solution.  If  a  simple 
MRF  is  used  to  model  the  data,  the  energy  function  will  yield  a  poor  disparity 
by  smoothing  across  object  boundaries,  particularly  when  occluding  objects 
are  present.  We  exploit  optical  flow  information  to  indicate  object  boundr^ies 
(depth  discontinuities)  and  occluded  regions,  in  order  to  improve  the  disparity 
solution  in  occluded  regions.  A  stochastic  relaxation  algorithm  (Simulated 
Annealing)  is  used  to  find  a  favorable  disparity  solution  by  minimizing  the 
energy  equation. 
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1.  INTRODUCTION 

Integration  of  visual  cues  have  recently  been  investigated  by  many  researchers  [1]-  [2]. 
Stereo  vision  [3]-[4]  and  optical  flow  [5J  information  have  been  used  to  obtain  3-dimensional 
structure  about  the  environment.  The  major  problems  encountered  in  stereo  marching  are 
identifying  the  occluded  areas  and  the  depth  discontinuities.  Our  objective  in  this  research  is 
to  develop  a  motion-stereo  matching  algorithm  that  is  minimally  affected  by  occlusion  and 
depth  discontinuities.  This  is  down  by  integrating  the  edge  information  and  short  range  opti¬ 
cal  flow  data  into  our  stereo  matching  algorithm. 

In  this  report  a  motion-stereo  vision  technique  is  investigated  where  short  range  optical 
flow  information  [5]  is  integrated  into  an  intensity-based  stereo  vision  technique.  A  Bayesian 
model  is  used  to  derive  the  Maximum  A  Posteriori  (MAP)  stereo  matched  solution  for  a 
motion-stereo  intensity  matching  algorithm.  The  input  data  and  the  disparity  solution  are 
modeled  as  Markov  Random  Fields  (MRFs).  The  MRF-Gibbs  Distribution  equivalence  [6] 
reduces  the  MAP  problem  to  that  of  finding  an  appropriate  energy  function  that  describes  the 
constraints  on  the  solution.  Stereo  algorithms  based  on  low-level  information  (intensity,  edges 
for  example)  have  difiiculty  with  occlusion  as  well  as  depth  discontinuities  because  they  have 
to  make  decisions  about  object  boundaries  in  order  to  find  the  correct  disparity  for  a  region 
that  is  not  visible  in  both  images.  If  an  object  boundary  is  not  pronounced,  the  disparity  solu¬ 
tion  for  the  occluded  region  will  usually  be  a  smoothly  changing  function  between  the  dispar¬ 
ity  of  the  occluded  object  and  the  visible  object.  In  this  report  the  optical  flow  data  is  used  to 
flag  the  potentially  occluded  regions,  and  when  a  region  is  identified  as  occluded,  a  discon¬ 
tinuity  in  the  disparity  solution  is  assumed.  However,  the  resulting  energy  function  is  no 
longer  convex  so  a  stochastic  relaxation  technique  is  employed  to  find  the  minima  of  the  func¬ 
tion.  In  stereo  images  occlusion  occurs  when  some  points  are  visible  in  one  image  but  not  in 
the  other  due  to  the  displacement  between  the  two  cameras.  For  instance,  this  will  happen 
when  an  object  looks  like  it  is  behind  another  one  from  the  angle  of  the  camera.  Figure  1 
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shows  the  mismatch  resulting  from  a  stereo  matching  algorithm  which  does  not  consider 
occlusion.  Occlusion  problem  has  been  tackled  by  other  researchers  who  have  used  the  rate  of 
change  of  disparity  [7],  or  the  "ordering  constraint"  [8]  to  flag  a  potentially  occluded  region. 

-In  Section  2  of  this  report  we  describe  our  algorithm  in  detail,  and  presents  properties  of 
MRPs  and  their  Gibbs  equivalence  distribution,  as  well  as  the  used  energy  functions.  In  sec¬ 
tion  3,  the  effect  of  occlusion  and  depth  discontinuities  are  reported  with  experimental  results. 

2.  DESCRIPTION  OF  THE  MOTION-STEREO  MATCHING  ALGORITHM 

The  approach  of  our  stereo  matching  algorithm  is  to  find  an  estimate  of  the  "optimal" 
disparity  solution  for  a  pair  of  stereo  images.  Optimality  is  defined  as  the  Maximum  A  Pos¬ 
teriori  probability  (MAP)  solution.  The  MRF  model  allows  us  to  integrate  the  low  level 
visual  cues  as  well  as  reducing  the  stereo  matching  problem  to  that  of  finding  an  appropriate 
energy  function.  We  minimize  the  energy  function  and  obtain  a  MAP  estimate  by  using  a 
well  established  stochastic  optimization  technique.  In  the  following  subsections  we  introduce 
MRFs  and  the  Coupled  MRFs,  MRF-Gibbs  distribution  equivalence,  energy  functions 
representing  the  constraints  on  the  solutions. 

2.1  Markov  and  Coupled  Markov  Random  Fields 

The  stereo  matching  problem  is  an  ill-posed  problem  in  the  sense  defined  by  Hadamaid 
[9-10],  which  can  be  solved  by  using  a  standard  regularization  method  [11-13].  The  standard 
regularization  theory  produces  a  convex  quadratic  energy  function  (a  cost  functional)  thus  with 
a  unique  solution.  However,  in  order  to  include  the  occlusion  process  and  to  deal  effectively 
with  the  integration  of  visual  cues,  a  non-convex  energy  function  is  derived.  Therefore,  in 
contrast  to  standard  regularization  methods,  we  represent  the  a  priori  knowledge  in  terms  of 
probability  distributions  rather  than  constraints  on  the  actual  solution  space.  Specifically,  we 
maximize  the  a  posteriori  probability  of  the  disparity  solution  given  the  degraded  (i.e.,  noisy) 
images  by  modeling  the  disparity  solution  as  a  Markov  Random  Field  D  over  a  lattice 
Assuming  that  in  general  the  disparity  is  smooth  over  rigid  objects  except  at  depth 
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discontinuities  the  MRF  property  can  be  stated  by  Equation  (1)  which  says  that  probability  of 
a  disparity  value  dij  at  a  pixel  location  (/,  j),  given  the  disparity  of  all  other  pixels  in  the  image, 
is  the  same  as  the  probability  of  that  value  given  only  the  disparity  of  the  neighbors  of  site 
O',  j) 

Probi  0=dij  I  da.  ikj)  ^  (ij))  =  Prob(  D  =  dij\  da.  (k,l)  e  N^iJ)  )  (1) 

where  NaiiJ)  specifies  the  pixel  indices  of  those  defined  to  belong  to  the  neighborhood  system 
of  o,y). 

The  disparity  data  with  the  depth  discontinuity  process  or  occlusion  process  can  be  com¬ 
bined  together  to  form  a  coupled  MRF  [6].  For  example,  define  the  depth  discontinuity  as  a 
binary  MRF  Ap,^0'jv.7iv)  defined  over  a  lattice  Us  as 

1 1,  if[p’’0'.7)-p''0w.7Ar)]^  ^  Tp 
~  I  0,  Otherwise 

where  is,  jstNpiij)  represents  the  eight  nearest-neighborhood  system,  Tp  is  a  pre-defined 
threshold,  and  p'ii,  j)  represents  the  magnitude  of  the  optical  flow.  This  binary  process  is  very 
similar  to  the  line  process  introduced  by  Geman  and  Geman  [6].  We  now  associate  this 
binary  process  with  the  disparity  process  to  form  a  coupled  MRF  defined  as  D  =  (D,  Ap)  over  a 
lattice  Hs  =  Vs\jUs  where  there  is  a  line  process  between  each  pair  of  disparity  sites. 

Each  input  image  can  also  be  considered  to  be  a  combination  of  two  MRF’s  coupled 
together  to  represent  the  corresponding  intensity  and  optical  flow  information.  For  example, 
let  the  intensity  image  be  a  MRF  X  on  a  lattice  Zs  such  that, 

X  =  /(«,  J)  if  ( iJ)  tZs  ;  Zw  =  I  (iJ) :  1  <  (iJ)  <  N  |  (3) 

where  /(/,  J)  is  the  gray  level  of  the  image  of  size  Nx  N. 

The  neighborhood  system  for  each  pixel  /(/,  j)  is  defined  as  its  eight  nearest-neighboring 
pixels.  Similarly  the  optical  flow  information  can  also  be  considered  as  a  MRF  O  on  a  lattice 
Ws  such  that. 
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O  =  p(i.  i)  if  ( i,j)  tW^  Wjv  =  I  iij) ;  1  <  OJ)  <  A/  I  (4) 

Where  p(i,  j)  represents  the  magnitude  of  the  optical  flow  at  location  (j,  j).  The  neighborhood 
system  Co  for  each  pixel  is  defined  as  the  eight  nearest-neighbors. 

Now  each  input  image  can  be  considered  as  a  coupled  MRF  F  =  (X,  O)  a  collection  of  an  inten¬ 
sity  process  X  and  optical  flow  process  O.  The  values  of  the  field  F  is  defined  as: 

=  if  iJeZs 

^  ^  \Oi,i  =p(i.  j)  if  ij  e  WV 

where  F  is  defined  as  a  MRF  on  a  lattice  the  components  of  F  assumes  values 

among  the  allowable  gray  level  and  optical  flow  values.  The  interaction  between  the  intensity 
process  X  and  the  optical  flow  process  O  will  become  clear  in  section  2.3  when  correspondence 
problem  (energy  expression)  is  considered.  We  can  define  a  pair  of  stereo  images  as  two  cou¬ 
pled  MRFs  f'  =  (X‘,  0‘)  and  F'=  (X^  O')  representing  the  left  and  the  right  image  data  respec¬ 
tively.  Each  coupled  MRF  can  be  considered  as  an  observation  or  an  external  input  to  the 
stereo  system. 

2.2  MRF-Gibbs  Distribution  Equivalence 

It  has  been  shown  [14-17]  that  if  an  image  measure  (such  as  disparity)  can  be  modeled 
by  a  MRF,  then  the  probability  distribution  of  the  measure  can  be  represented  as  a  Gibbs  dis¬ 
tribution  (Hammersley-Cliffbrd  theorem).  Equation  (6)  defines  the  probability  of  arriving  at  a 
solution  configuration  (o,  with  a  Gibbs  distribution  given  by 


P(D=(o)  = 


,-Uia>)/kT 


where  Z,  is  the  normalizing  constant  such  that  the  sum  of  P{(x))  for  all  configurations  is  unity: 

Z  =  '£cxpl-U((o)/kT].  ( 


T  is  analogous  to  the  temperature  of  the  system  that  yields  the  configuration  co,  and  i  is  a 
universal  constant.  C/(o))  is  called  the  energy  function  and  is  of  form. 
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f/((o)=5:  f/,((o)  (8) 

c 

where  U^iai)  is  called  a  potential  function  and  is  defined  over  the  neighborhood’s  cliques.  A 
clique  is  a  set  of  sites  (pixels)  such  that  for  a  defined  neighborhood  system  every  pair  of  sites 
belonging  to  this  set  are  neighbors.  An  example  of  neighborhood  system  and  possible  poten¬ 
tial  functions  used  in  this  proposal  is  given  below: 

1)  The  potential  function  for  a  simple  disparity  MRF  with  no  discontinuity  process  at  a 
pixel  location  (/,»  with  a  neighborhood  system  is  defined  by 

VMj)  =  £  [  j)  -  d  Us,  jN)f  isJN  e  m,  j)  (9a) 

‘mJh 

where  the  neighborhood  system  Na(i,  j)  represents  only  the  eight  nearest-neighbors. 

2)  The  potential  function  for  a  coupled  MRF  (disparity  with  discontinuity)  at  pixel  location 
(i,j)  with  a  neighborhood  system  j)  is  defined  by 

^Mj)  =  X  t  -  diiuJM)  f  •  (l-ApiyOW Jm))  +  •  Cp;  imJm  e  Mdih  j)  (9b) 

in -in 

where  j)  is  simply  the  8-nearest  neighbors  and  Cp  is  a  constant.  Equation  (6),  with  the 
potential  function  given  by  expression  (9a)  or  (9b)  is  said  to  represent  the  prior  distribution  of 
the  disparity  solution.  It  is  easily  seen  that  low  energy  configurations  represent  a  higher  pro¬ 
bability  than  high  energy  configurations.  The  ^=kT  factor  affects  the  equilibrium  state  of  the 
solution.  At  high  temperatures,  the  solution’s  equilibrium  state  will  be  very  random.  At  low 
T,  however,  the  equilibrium  state  is  more  regular  and  dependent  on  its  initial  configuration. 

The  advantage  of  representing  each  probability  in  the  Gibbs  distribution  form  is  that  it 
can  1 formulated  with  an  energy  function  U  (oa)  which  suitably  describes  the  constraints  of  the 
solution.  These  constraints  are  similar  (if  not  identical)  to  those  used  in  a  standard  regulariza¬ 
tion  method,  but  the  interactions  between  the  depth  discontinuity  and  disparity  smoothness 
processes  can  now  be  explicitly  included  in  the  energy  function. 

2.3  The  Energy  Function  for  the  Proposed  Motion-Stereo  Matching  Algorithm 
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The  energy  function  describes  the  constraints  on  the  desired  solution  in  terms  of  local 
characteristics  [18].  For  example,  there  should  be  a  strong  similarity  in  the  intensity  (or  orien¬ 
tation  if  an  edge)  between  two  corresponding  sites  in  the  stereo  solution.  Likewise,  the 
disparity  solution  should  be  smooth  over  the  same  object. 

The  MRF  model  allows  us  to  reduce  the  stereo  matching  problem  to  that  of  finding  an 
expression  for  the  a  posteriori  probability  P(D  IF',  F^)  in  terms  of  the  a  priori  probability  F(D) 
and  the  conditional  probability  of  observations.  Then  using  the  MRF-Gibbs  (or  Boltzmann) 
distribution  equivalence,  one  can  determine  an  appropriate  energy  function  representing  the 
corresponding  MAP. 

The  left  and  the  right  image  data  are  defined  as  two  observations  corresponding  to  ran¬ 
dom  fields  f'  and  F'  respectively.  Our  goal  is  to  determine  the  most  likely  solution  D  =  (d,_y) 
given  the  above  observations.  Using  the  Bayesian  rule  the  a  posteriori  probability  can  be  writ¬ 
ten  as 


PiD\F'  F(F',F-.£>)  P(F'|£>,F0-F(£>,F0  P{F'\D,  F')  ■  P{D\r) 

P(F‘,F'-)  P{F\F')  F(F'IF0 

If  the  disparity  D  is  independent  of  the  observation  F'  then  P(p\F')  =  P{D)  which  can  be 
replaced  by  its  corresponding  Gibbs  distribution  given  by  expression  (6).  The  conditional  pro¬ 
bability  F(F'  I  F')  is  a  constant  since  it  is  independent  of  D.  The  key  term  in  expression  (10)  is 
the  conditional  probability  F(F'lD,  F*")  which  has  to  be  written  in  Gibbs  distribution  form. 
From  the  stereo  camera  geometry  we  have  F'{i,  j)  =  F'(i,  j+dij)  in  absence  of  any  noise  or 
occlusion. 


However,  assuming  F'  was  degraded  by  a  Gaussian  random  noise  N  independent  of  D  and 
F'  with  zero  mean  then. 


20^ 


(11) 


1 


p{f'\d,f')=  n  , — 


(1%  j+d,i)f 


2(P- 


(12) 
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Expression  (12)  can  be  written  in  Gibbs  distribution  form  as 

P(F‘\D,F')=^cxp^^-^  £  t/i/F'lD.F")  J  (13) 

where 

Uij(F‘  ID,  F')  =  ^  [( /'(/.  j)  -  m,  j+d,j)  )2  ]  (14) 

with  X.= 

20^ 

The  overall  local  energy  contribution  from  a  single  site  (i,  j)  can  be  written  as 

=  t/y(F'  ID.  FG  +  VMij)  (15) 

where  the  first  term  represents  data  term  the  similarity  in  the  intensity  (or  orientation  if  an 
edge)  between  two  corresponding  sites  in  the  stereo  solution  as  given  by  expression  (14).  The 
second  term  represents  the  smoothness  of  the  solution  given  by  expression  (9a)  or  (9b)  when 
depth  discontinuity  is  present.  Tbe  total  energy  is  determined  by  tof’iling  the  energy  contri¬ 
bution  from  each  site. 

In  order  to  include  depth  discontinuities  as  weU  as  occlusion  process  we  assume  that  for 
each  site  the  optical  flow  estimate  is  available  for  the  left  and  the  right  images.  There 
are  a  number  of  standard  techniques  that  can  be  used  to  measure  the  optical  flow  information 
[5].  Optical  flow  is  represented  by  a  vector  random  field  where  ^/,y)  =  («;; .  Vy),  the  two  com¬ 
ponents  of  the  optical  flow.  The  difference  in  optical  flow  between  the  two  corresponding 
sites  is  thresholded  to  indicate  a  potentially  occluded  region  (4)p  =  l)  or  a  visible  region 
(<l>p  =  0).  only  the  magnitude  of  the  optical  flow  is  considered. 

J  I,  iff  lp'(/,y)l  -  \p'(i+dij,j)\ 

~  0,  Otherwise 

An  occluded  region  is  then  processed  differently  than  visible  regions;  heavily  weighting 
the  smoothness  term  over  the  data  term.  As  shown  in  Equation  (17),  the  data  term  Dy  does 
not  use  the  measure  of  similarity  of  the  corresponding  pixels  in  an  occluded  region  as  they 
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will  not  match  for  the  correct  disparity  value.  Instead,  the  data  term  is  assigned  a  small  fixed 
energy  causing  the  smoothness  term  to  have  more  weight  in  determining  the  correct  dispar¬ 
ity  estimate.  The  non-zero  cost  (energy)  associated  with  this  process  inhibits  non-occluded 
areas  firom  settling  into  a  falsely-indicated  occluded  state. 

where  (j,))  e  Sn-  The  smoothness  term  Si,  is  not  changed  significantly  except  that  we  now  use 
the  optical  flow  gradient  to  indicate  potential  discontinuities  in  the  disparity  solution.  If  the 
difference  in  optical  flow  of  two  neighboring  sites  exceeds  a  pre-specified  threshold  Tp  then 
there  is  a  discontinuity  in  the  disparity  and  smoothness  measure  should  not  be  performed 
across  this  discontinuity.  Instead  a  small  energy  cost  Cp  is  included,  in  this  way  we  are  intro¬ 
ducing  a  local  minima  in  our  cost  function. 

Sij  =  2  f  ^0«y)  ~  ^ {iuJu)  ■  (l-APv(*w./i/))  +  ■  Cp;  iuju  t  MjiiJ)  (18) 

‘mJm 


where 


•  1,  if  [  lpTiJ)i  -  1pT/mJm)I  ]"  S  T, 
0,  otherwise 


(19) 


If  the  calculated  optical  flow  is  reliable  it  can  be  segmented  into  different  regions  representing 
each  object  and  their  boundaries  [19]  or  a  cooperative  technique  could  be  designed  that  can  do 
the  segmentation  as  the  optical  flow  is  calculated. 


3.  EXPERIMENTAL  RESULTS 

We  work  with  motion-stereo  images  generated  by  two  cameras  situated  in  a  parallel  axis 
geometry  yielding  epipolar  lines  which  are  parallel  to  the  scan  lines.  Thus,  corresponding 
image  points  lie  along  the  same  horizontal  scan  line  limiting  the  search  space  to  one  dimen¬ 
sion.  For  the  initial  algorithm  validation,  we  simulated  the  intensity  and  optical  flow  vector 
for  stereo  images  of  size  64  x  64  pixels.  Our  simulated  image  model  includes  random  errors 
in  both  the  image  intensity  and  the  optical  flow  estimate.  Figure  2  shows  a  pair  of  stereo 
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images  simulated  for  experimental  evaluation  of  the  proposed  technique. 

For  the  preliminary  experimental  investigation  it  is  assumed  that  optical-flow  values  for 
the  motion  stereo  images  are  already  estimated  from  a  motion  estimation  technique.  Although 
there  are  problems  with  calculating  reliable  optical  flow  data,  but  by  means  of  post-processing, 
the  optical  flow  information  can  be  smoothed  and  segmented  into  regions  representing 
different  objects  and  can  be  used  to  locate  the  occluded  regions  as  discussed  in  this  report. 

To  maximize  the  probability  given  by  Equation  (10),  and  thus  obtain  the  a  MAP  esti¬ 
mate,  we  minimize  the  total  energy  by  achieving  an  equilibrium  state  at  low  temperature.  We 
use  a  stochastic  relaxation  method  called  Simulated  Annealing  (SA)  which  is  described  in 
detail  in  [20-21].  The  physical  process  of  annealing  brings  materials  to  low  energy  states  by 
gradually  lowering  their  temperature.  As  in  chemical  annealing,  if  the  temperature  is  dropped 
too  quickly  the  solution  does  not  have  sufficient  time  to  reach  equilibrium  states,  resulting  in 
configurations  which  correspond  to  local  minima  of  the  energy  function;  these  constitute 
errors.  The  ideal  temperature  schedule  which  guarantees  convergence  to  a  minimum  energy 
state,  or  the  maximum  a  posteriori  (MAP)  configuration  is  given  in  [6]  if  the  Gibbs  sampler 
schedule  is  used;  however,  this  schedule  is  impractically  slow.  For  the  preliminary  evaluation 
of  the  algorithm  we  have  used  a  non-ideal,  monotonically  decreasing  temperature  schedule 
and  some  result  are  presented.  The  Simulated  Annealing  algorithm  can  efficiently  be  imple¬ 
mented  with  parallel  processing  units,  such  as  neural  networks,  due  to  its  localized  nature. 
The  stochastic  characteristic  of  the  SA  algorithm  however,  requires  a  random  parameter  to 
simulate  the  effect  of  lowering  the  temperature. 

Figure  3  illustrates  the  disparity  solution  when  no  occlusion  process  is  employed.  The 
energy  function  used  is  a  combination  of  the  solution  smoothness  and  intensity  similarity  con¬ 
straints  given  by  expressions  (9a)  and  (14)  respectively.  Disparity  error  occurs  in  the 
occluded  areas  when  this  energy  function  is  used.  Here  the  solution  has  matched  to  the  wrong 
pixel  in  an  attempt  to  minimize  the  data  term.  Figure  4  shows  the  result  of  applying  Equa¬ 
tions  (17)-(18)  to  images  similar  to  those  in  Figure  2.  The  solution  is  much  improved  and 
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yields  a  better  disparity  estimate  in  the  excluded  region.  Figure  5  shows  disparity  solution  at 
three  different  temperatures,  T  =  50,  10,  and  l  for  the  intensity-based  stereo  technique  incor¬ 
porating  optical  flow  data.  The  temperature  was  decreased  by  10%  when  there  was  no 
significant  change  in  the  energy  value  (system  in  an  equilibrium  state).  The  system  is  said  to 
be  in  an  equilibrium  state  whenever  the  number  of  jumps  in  the  energy  increase  and  decrease 
are  about  the  same.  Although  there  is  no  guarantee  of  convergence  to  a  minimum  energy  state 
due  to  our  monotonic  temperature  schedule,  but  a  close  approximation  to  the  optimal  solution 
is  expected  whenever  the  system  is  kept  in  equilibrium  at  each  temperature. 
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Figure  1.  Typical  stereo  matching  error  in  occluded  regions. 


Figure  2.  Typical  left  and  right  stereo  images  used  to  generate  results  in  Figures  3-7,  Gaus¬ 
sian  noise  is  also  introduced  into  the  data.  Disparity  of  left  box  =15;  of  right  box  = 


